
function [Wf,Wi,Wn,iterVF]=solve_value_function_A1(Par1,Ff_1,wf_1,Fi_1,wi_1,n)

r=0.036;

% Clenshaw-Curtis quadrature
% n+1 nodes: x in [-1,1] by decreasing order
[x,wcc] = ccquad(n-1); x=flipud(x);

transition_rates=Par1(1:6)
theta=Par1(11)
a=Par1(12)
gamma=0; % gamma = 0 at benchmark
b1 = Par1(13)
df_1=transition_rates(1);
di_1=transition_rates(2);
lnf_1=transition_rates(3);
lni_1=transition_rates(4);
lfi_1=transition_rates(5);
lif_1=transition_rates(6);


%preallocation to speed up

Wi = zeros(n,1); Rw_f_i = zeros(1,n); Rw_i_f = zeros(1,n); Wf_Rw_f_i = zeros(1,n); Wi_Rw_i_f = zeros(1,n);
Wf = zeros(n,1); Fi_1_Rw_i_f = zeros(1,n); Ff_1_Rw_f_i = zeros (1,n); derWf = zeros(1,n);derWi = zeros(1,n);
intmaxWf_Wi = zeros(1,n); intmaxWi_Wf = zeros(1,n); eWi = zeros(1,n); eWf = zeros(1,n);

%% Transforming the probabilities g and f in densities
% If wages are linearly spaced
wht_f_1=(wf_1(n)-wf_1(1))/(n-1); 
wht_i_1=(wi_1(n)-wi_1(1))/(n-1);


    % Initializing
    
    % value function (min and max)
    Wn=(1-exp(-theta*(min(wi_1))))/theta/r;
    Wf_min= (1-exp(-theta*(min(wf_1))))/theta/r; Wf_max= (1-exp(-theta*(max(wf_1))))/theta/r;
    Wi_min=(1-exp(-theta*(min(wi_1))))/theta/r; Wi_max=(1-exp(-theta*(max(wi_1))))/theta/r ;
    
     
    Wf=(Wf_min+Wf_max)/2+x*(Wf_max-Wf_min)/2; % Interpolation with CC quadrature (defined above)
    Wi=(Wi_min+Wi_max)/2+x*(Wi_max-Wi_min)/2;

    
    critVF=1;
    iterVF=1;
    omegaW=0.1;
    while mean(critVF)>0.005 && iterVF<=400
        
        % Reservation wages
        for i=1:n
            % to construct Wfn
          
            dWi_Wf=abs(Wi-Wf(i))./Wf(i);
            Rw_i_f(i)=wi_1(min(find(dWi_Wf<=min(dWi_Wf))));
            Wi_Rw_i_f(i)=Wi(min(find(dWi_Wf<=min(dWi_Wf))));
            Fi_1_Rw_i_f(i)=Fi_1(min(find(dWi_Wf<=min(dWi_Wf))));
        
            
            derWf(i)=1./(r+df_1 +lfi_1*Fi_1_Rw_i_f(i));
            
            % to construct Win
            dWf_Wi=abs(Wf(i,:)-Wi(i))./Wi(i);
            Rw_f_i(i)=wf_1(find(dWf_Wi<=min(dWf_Wi), 1 ));
            Wf_Rw_f_i(i)=Wf(find(dWf_Wi<=min(dWf_Wi), 1 ));
            Ff_1_Rw_f_i(i)=Ff_1(find(dWf_Wi<=min(dWf_Wi), 1 ));
            
            
            derWi(i)=1./(r+di_1+lif_1*Ff_1_Rw_f_i(i));
        end
       
                % to construct Wnn
                dWf_Wn=abs(Wf-Wn)./Wn;
                Rw_f_n=wf_1(find(dWf_Wn<=min(dWf_Wn), 1 ));
                Ff_1_Rw_f_n=Ff_1(min(find(dWf_Wn<=min(dWf_Wn))));
                
                dWi_Wn=abs(Wi-Wn)./Wn;
                Rw_i_n=wi_1(find(dWi_Wn<=min(dWi_Wn), 1 ));
                Fi_1_Rw_i_n =Fi_1(min(find(dWi_Wn<=min(dWi_Wn))));
        
        % solve integral by parts
        
        % Wf
        for i=1:n
            intmaxWi_Wf(i)=sum((wi_1>Rw_i_f(i)).*Fi_1.*derWi'.*wht_i_1); 
        end
        % Wi
        for i=1:n
            intmaxWf_Wi(i)=sum((wf_1>Rw_f_i(i)).*Ff_1.*derWf'.*wht_f_1);   
        end
        
        % Wnn
        
        intmaxWf_Wn=sum((wf_1>Rw_f_n).*Ff_1.*derWf'.*wht_f_1);
        intmaxWi_Wn=sum((wi_1>Rw_i_n).*Fi_1.*derWi'.*wht_i_1);
        
        %%%%% Value functions
        % solve in terms of F, reservation wages and lambdas
        
        % Wfn(w1) and Win(w1)
        for i=1:n
            
            eWf(i)=(1/(r+df_1))*((1-exp(-theta*(wf_1(i))))/theta+a+df_1*Wn ...
                +lfi_1*intmaxWi_Wf(i))-Wf(i);
            
            eWi(i)=(1/(r+di_1))*((1-exp(-theta*(wi_1(i))))/theta+gamma+di_1*Wn ...
                +lif_1*intmaxWf_Wi(i))- Wi(i);

            
        end

        % Wnn
        
        eWn=(1/r)*((1-exp(-theta*b1)))/theta+gamma ...
            +lnf_1*intmaxWf_Wn ...
            +lni_1*intmaxWi_Wn;

        
        % Update value functions
        
        Wf=Wf+omegaW*eWf';
        Wi=Wi+omegaW*eWi';
        Wn=Wi(1);eWn=0;
        
        critVF=[abs(eWf'./Wf);abs(eWi'./Wi);abs(eWn/Wn)];
        
        iterVF=iterVF+1;
     
    end
    iterVF

    Rw_1_f_ind = Rw_f_n; 
    Rw_1_i_ind = Rw_i_n;  
    save Rw_1_f_ind Rw_1_f_ind
    save Rw_1_i_ind Rw_1_i_ind
end
 